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TECHNICAL PUBLICATION 


A RECOMMENDED PROCEDURE FOR ESTIMATING THE COSMIC-RAY SPECTRAL 
PARAMETER OF A SIMPLE POWER LAW WITH APPLICATIONS 

TO DETECTOR DESIGN 


1. INTRODUCTION 


This Technical Publication (TP) develops and compares two statistical methods for estimating 
the spectral parameter of the simple power law energy spectrum from simulated detector responses 
(energy deposits). The maximum likelihood (ML) procedure, which is shown to be the superior 
approach, is then generalized for application to a set of real cosmic-ray data that make the methodol- 
ogy applicable to existing cosmic-ray data sets. 

As part of this research, analytical methods were developed in conjunction with a Monte Carlo 
simulation to explore the combination of the expected cosmic-ray environment with a generic 
space-based detector and its planned life cycle. This allows exploration of various detector features 
and their subsequent impact on estimating this spectral parameter. This study thereby permits instru- 
ment developers to make important trade studies in design parameters as a function of the science 
objectives, which is particularly important for space-based detectors where physical parameters, such 
as dimension and weight, impose rigorous practical limits to the design envelope. 


2. SIMPLE POWER LAW 


The simple power law suggests that the number of protons detected above an energy E for an 
assumed collecting power (combination of size and observing time) is given by 1 


N 0 OE) = N a 


< £ V a i +I 
\ E AJ 


0 ) 


where E is in units TeV, a | is believed to be ~2.8, and N A and E A are numbers determined from the 
detector size and exposure time in the environment. For a typical space-based detector of 1 m 2 with a 
3-yr program life, N A and E A are 160 and 500 TeV, respectively, implying that this detector is 
expected to observe 160 proton events above 500 TeV over its expected life cycle. In statistical terms, 
N 0 is assumed to represent an average number of events while the actual number to be observed on 
any given mission would follow the Poisson probability distribution with mean number N () . The num- 
ber of particles detected depends only on the geometrical factor of the assumed detector and its mate- 
rial composition. The detection efficiency is a convolution of the geometry and material composition 
and is taken to be independent of energy. 

The associated cumulative probability distribution function for£ over an energy range [£j^ 2 1 
is then given by 


<V £ > 


N q (> E) — N q (> £) ) 
N 0 (> E^) — Nq(> Ej ) 


for fj < E < E 0 


= 1 


E ~ a+ 1 — e\ a ' 

£ ,_a 1 _ £* _0: 1 
c, c 2 


(2) 


Thus, the corresponding probability density function (pdf) for £ is obtained by differentiating 
equation (2) to give 


4 > o( £ ) = 


^p(E) 

dE 
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for E\< E < E-, 


(3) 
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To randomly sample GCR proton event energies from the simple power spectrum over the 
interval [E X ,E 2 ], w~4> 0 (£,) is solved in terms of £, to obtain 


1 

E i =% 1 '(«,■) = 

where </, is a simulated random number from a standard uniform distribution and <J> 0 _1 represents the 
inverse function of <D 0 , which is a conventional notation that will be used in subsequent sections. The 
mean, variance, and other moments of the simple power law distribution are determined by the expected 
value operator, where the general form of (£"') is 
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Note a crucial point at this time, and that is <E 2 > becomes infinite (as do all other higher 
moments) as £7 goes to infinity, which is easily seen in equation (6): 


lim x 2 x ^ dx = 00 for all A < 3 and a > 0 . ^ 

b-> oo J <> 

This observation suggests the need for a careful look at the effects of the large variance and 
other higher moments associated with all power law distributions, even when E 2 is kept finite. A 
measure of the relative dispersion of the energies of the incident protons, which is independent of 
units, is defined by V=o E /p E for the simple power law and is called the coefficient of variation in the 
statistical literature. An important concept in detector design is the energy resolution p of the detector 
that provides a measure of the relative accuracy of a cosmic-ray detector, which is the fractional error 
in measurements of a monoenergetic beam. The resolution p is defined as the standard deviation 
divided by the mean response with typical values of 30 to 40 percent. 

As will be shown in this TP, the precision with which the spectral parameter a, can be esti- 
mated from a set of detector responses (energy deposits), measured in terms of its standard deviation, 
is a function of both the variance of the incident energies and the uncertainty induced by the detector. 
The dominating component of this measurement precision will be shown to be attributable to the 
standard deviation of the incident energies o E , which in turn can only be controlled through collecting 
power. Since V and p are dimensionless and provide a measure of relative dispersion for the power 
law distribution and detector, respectively, an instructive comparison will show that V » p. To illus- 
trate these points, a detector life cycle having parameters N A = 160 and E A = 500 TeV will observe 
52,200 events on average in the energy range £, = 20 TeV to £ 2 = 5,500 TeV from a simple power law 
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spectrum when a, is 2.8. This gives a mean galactic cosmic-ray (GCR) event energy p E = 44.5 TeV. 
a standard deviation a E = 74. 1 0 TeV, and a coefficient of variation V = 1 66.5 percent. In comparison, 
the resolution p of most detectors is between 30 and 40 percent. E-, is chosen for this detector life 
cycle combination as 5.500 TeV, since the expected number of events above this energy are negli- 
gible, while E | is taken to be 20 TeV for purposes of this discussion. 

Since the number of events and their incident energies will vary because of the finite detector 
size and exposure time, the statistical behavior of the GCR event energies in combination with a 
detector having energy resolution p and the subsequent spectral parameter estimate over multiple 
missions shall be studied. Thus, for each mission, a random number N of GCR events from a Poisson 
distribution with mean 52,200 representing the number of simulated events that the detector will 
observe in the energy range 20 to 5,500 TeV on any given mission is first generated. 

Next, the incident energy of each of these N events using equation (4) is simulated. For 
example, for one such simulated mission, N = 51,883 and the mean and standard deviation of the 
simulated GCR incident energies are calculated to be 43.85 and 66.39 TeV, respectively. To illustrate 
the large fluctuations associated with power law distributions, the same number of events (51,883) 
from a normal distribution having a mean of 44.5 and standard deviation 74.1 so as to match the 
power law’s mean and standard deviation for this energy range when a j = 2.8 was also simulated. The 
sample mean and standard deviation, 44.5 1 and 74. 1 7, respectively, for a single sample mission, which 
are much closer to the population mean and variance than those from the power law random samples, 
were also observed. The process is repeated for 100 missions and the standard deviation for each 
mission is plotted in figure 1. 


- Power Law Sigma 

- Sigma (Normal Distribution) 



Figure 1 . Standard deviation of simulated incident energies from power 
law (ragged curve) for 100 missions compared with that from 
normal distribution having same mean and variance. 
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Note the large fluctuations of the standard deviations for the power law samples from mission 
to mission; while in contrast, the standard deviations of missions generated from a normal distribu- 
tion are very stable. As will be seen in subsequent sections, this is why the variation in detector 
responses is dominated by the variation of the GCR event energies, relative to the uncertainty induced 
by detector resolution. This in turn contributes the dominant component of the standard deviation of 
the spectral parameter estimator. 

The variation of the sample standard deviation ,v used as an estimator of a and measured by its 
standard deviation is given by 


= (7) 

* ] 4jU 2 N ’ 

where /t,. is the rth central moment about the mean, 2 defined for the simple power law as 

IU r =f(E-u, : Y(l) 0 (E)dE . (8) 

Thus, the large variation in mission standard deviations is due to the term fJ 4 , which again is 
only finite by setting to a finite value, but nevertheless is responsible for the erratic behavior of the 
mission-to-mission sample standard deviations depicted in figure 1. This erratic behavior of the 
observed mission standard deviations will be true for any power law having spectral index a { in the 
range 3 <a\ <5. Note that for the normal distribution. 


a 


a 


S 


V2/V ’ 


(9) 


and evaluation of these two formulae yield a s = 5 TeV for the simple power law and 0.229 TeV for the 
normal distribution, which is roughly a factor of 22. 
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3. ESTIMATION OF THE SPECTRAL PARAMETER a, 


Of particular interest in the study of cosmic rays is the estimation of the spectral parameter a j 
from a set of data. In this section, two statistical procedures for estimating a j are developed and 
compared. Even though in practice the actual incident particle energies are never observed (but only 
a measure of their energy deposition from their passage through the detector), it is important to con- 
sider the concept of an ideal detector having zero resolution. Thus, such a detector would measure the 
GCR event energies exactly. 


3.1 Method of Moments 

The method of moments consists of equating the sample moments with the population 
moments, which in general leads to k simultaneous nonlinear algebraic equations in the k unknown 
population parameters. For the simple power law, there is only one parameter to be estimated, so the 
sample mean£ is set to the population mean n E obtained by setting m - 1 in equation (5) to obtain 
equation ( 10) to be solved in terms of a, : 


£ = 


{ a, - 1 

V“i ~ 2 ) 


-a, 


- £ 


2-a, 




GO) 


Thus, for a given sample of size N, this equation is solved in terms of ttjby numerical methods 
to provide an estimate of a t . This estimator, which is a function of the random variable E . has its own 
associated pdf. Since the GCR incident energy E has mean fi E and finite variance o\ (only because 
the upper energy £ 2 is finite), it is known by the central limit theorem that the distribution of the 
sample average E follows a normal distribution with mean fi E and variance aim. 

For example, when oq = 2.8, £| = 20 TeV, and £ 2 = 5,500 TeV, £ is normally distributed with 
mean 44.5 TeV and standard deviation (74.1 T eV)/(V ,/2 . These results can be used to obtain the 
probability distribution of the estimator by solving the probability equation, 


Pr 


( a ] - 1 


-a, 


- E 


2-a, 


e\ a ' - e\ 


- a , 


74.1 

7 W 


44.5 

< Z > 



—oo 


2 d.X , 


GO 


6 


in terms of a, for various values of Z. Letting Z vary from -4.7 to 4.7 and setting N - 52,000 events 
gives the probability distribution of a, shown in figure 2. Also depicted in figure 2 is the relative 
frequency histogram of the estimates ct| , based on 5,000 simulated missions. For each mission, 52,000 
events on average are simulated and the estimate of a, obtained by solving equation (11). Further- 
more, even though there is no explicit mathematical form for the pdf, its mean and standard deviation 
can be calculated by numerical methods. For the distribution shown here, its mean is numerically 
evaluated to be 2.800 and standard deviation as 0.01 1 5 when N = 52,000, which compares to the mean 
and standard deviation of the 5,000 simulated estimates with 2.800 and 0.01 14, respectively. With the 
ability to numerically construct this estimator’s pdf and moments, the important result found was that 
its variance is inversely proportional to the sample size N , which happens to be true for many common 
estimators (e.g., the sample mean, standard deviation, and median). For example, if the number of 
events is doubled, then the variance is halved; and it the number of events is halved, then the variance 
doubles. Note that this relationship between sample size and the standard deviation of the estimator 
a, is based on keeping £, and E 2 fixed, so that in practice, increasing the collecting power can reduce 

the variance. 


Theoretical Distribution of 3 Versus 
Histogram of Estimates Based on 5,000 Missions 



Figure 2. Probability distribution of method of moments estimate 

of a, with relative frequency histogram of spectral parameter 
estimates obtained from simulation. 
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3.2 Method of Maximum Likelihood 


The likelihood function of a random sample from the simple power law. regarded as a function 
of the single unknown parameter a j, is 


L(a,) = 


a | - 1 


-a. 


-E' 


\ 

N 

N 3 

-a, 


m 

) 


u=i y 


-a. 


^\ — Ej — E 2 


( 12 ) 


The method of ML seeks as the estimate of a x that value (say, a ML ) which maximizes the 
likelihood function so that L(a ML )> L(a, ) for all a x . Statistically, this means that the ML estimator 
leads us to a choice of a | that maximizes the probability of obtaining the observed data. In practice, it 
is often simpler to work with the logarithm of the likelihood function and seek solutions of (log L)' 
= 0 for which (log L)" < 0 (indicating a maximum), where the prime and double prime indicate the 
first and second derivative, respectively. Thus, equation (13) is numerically solved in terms of a x to 
obtain the ML estimate a ML : 


3logL 

3a, 


N 

a x - 1 


N 


(log £,)£,' -(log£ 2 )ri 


1-a, 


N 

5>g£-,=0 . 


/=l 


(13) 


The second derivative of the log-likelihood function is next obtained. Note that (log L)"< 0 for 
all a,, indicating that log L is concave; hence, there is a unique maximum, which was graphically 
observed by plotting log L as a function of a x : 


3 2 logL 
da} 


= —N 


(a, -l ) 2 


+ 


' E^E\ +a ^(\o g E 2 -\o % E x f 

\ 

{E 2 E^ - E x E 2 1 ) 2 

/ 


(14) 


By the Cramer-Rao inequality, the lower bound of the variance of any estimator a of a, is 
given by 


Var(a)>— ^ , (15) 

3 2 logL 

da} 

which is asymptotically attained by the ML estimator. Also note that it is inversely proportional to the 
number of events N as was the variance of the estimator obtained using the method of moments. Other 
important properties of ML estimators are that they are (I) asymptotically normally distributed and 
(2) consistency or asymptotically unbiased. Thus, a key question is, “For what values of N are these 
asymptotic properties achieved by the ML procedure?” 
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Based on the same 5,000 mission set discussed in section 3, the mean and standard deviation 
of the 5,000 ML estimates are 2.800 and 0.00782, respectively. Using equations (15) and (16), the 
Cramer-Rao bound is computed to be 0.00786 when N = 52,000 and a ML = 2.800, which compares 
very well with the simulation results. Furthermore, the frequency histogram of these 5,000 ML esti- 
mates resembled the normal distribution as stated in ( 1 ) of the above paragraph. A separate simulation 
study was conducted in which the sample size N was gradually reduced from 52,000 to 200 and the 
two asymptotic properties — attaining the Cramer-Rao bound and consistency were achieved b\ the 
ML estimates until N equals =2,000. A bias on the high side of an ^ f a '' ure t0 attain the Cramer- 
Rao bound became more and more evident as the number of events N diminished from 2,000 to 200. 


Another very important comparison is the ratio of the standard deviation of a ML to that of the 
estimator, obtained using the method of moments. Direct calculation shows this ratio is roughly 1 .45, 
implying that the ML procedure is significantly better than the method of moments when dealing with 
the simple power law. This result is not too surprising, however, because ML estimators in general 
have better statistical properties than the estimators obtained by the method ot moments.' 
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4. DETECTOR RESPONSE FUNCTION 


Based on GEANT simulations of energy deposition for monoenergetic protons at specified 
energies of 0.1, 1, 10. 100, 1,000, and 5,000 TeV, the Gaussian distribution provided a reasonable 
description of the distribution of energy depositions at each of these incident energies. 4 Furthermore, 
the mean detector response was well approximated by a linear function of incident energy in the range 
of interest for this study, typically between 10 and 5,500 TeV. Other detector response functions, such 
as a gamma distribution and another response function constructed from a combination of normal 
distributions having different parameters, will also be presented. 

The random variable Y is introduced to represent the detector’s response in terms of energy 
deposition of a GCR proton of incident energy E. The conditional mean response and standard devia- 
tion of Y for a particular incident energy E are modeled as /dy^ = {a + bE) and <7y\£ =(c + dE), 
respectively, where the four coefficients a , b, c, and d are estimated using linear regression in the 
GEANT simulation results. Thus, for each simulated incident GCR proton energy £■, the detector 
response is simulated as 


or 


Y i ~ E\\E. + a Y\E j Z i 


(16) 


Y i ={,a + bE j ) + (c + dE j )Z j , ( |7 ) 

with the nonnegativity constraint Y f > 0 and where Z, is a standard normal random number having zero 
mean and unit standard deviation. Thus, the detector response function is defined as 


g{y I E) = 



O' ^vl/: * 


V >0 


(18) 


where tj vj £ is a normalizing coefficient related to the truncation of the normal distribution resulting 
from the constraint y > 0. It is worth noting for constant resolution studies in which a Gaussian 
response function is assumed and p = cr/p is set to values 0.4 and 0.6, the corresponding detector 
energy resolution is 39 and 5 1 percent, respectively, and is rounded to 40 and 50 percent in the figures 
and tables in this TR 
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Thus, T] v | f is determined from 





Py\E 


( 19 ) 


where the lower limit of integration is -1 divided by the resolution function, given as 

p v |£ = ®y\E / Py\E = + dE) / (u + bE) . 

First, it is worthwhile to consider a detector having energy resolution p yAE = <J YiE fry ie 
a constant pand independent of the cosmic-ray’s energy E so that o Y \ E ~ P Py\E ’ w ^ ere P' ca ^ va ' ues 
of interest for p are 0, 0.2, 0.3, 0.4. and 0.6. It should also be noted that the normalizing coefficient r\ 
in equation (19) is constant whenever the detector resolution p is energy independent. 

Second, a case where fl Y E and O y , e are linear but their ratio is not a constant, so that the 
detector’s resolution is a nonlinear function of incident energy E, was investigated. For this second 
scenario, two studies were conducted in which the resolution is getting better, from 40-percent 
resolution at 20 TeV to 30-percent resolution at 5,500 TeV, and then getting worse, from 30-percent 
resolution at 20 TeV to 40-percent resolution at 5,500 TeV. These two energy-dependent cases are 
presented in section 8. 



5. PROBABILITY DISTRIBUTION OF THE DETECTOR RESPONSE 


The probability distribution for the detector response in the presence of the simple power law 
energy spectrum over the energy range [£|,£ 2 1 is 


E 2 

£ 0 (y;a|)= J g(y\E;p)<f> 0 (E;a l )dE , v>0 . 

£ i 

The spectral parameter a | has been explicitly included in the argument list of both the simple 
power law pdf as <pQ(E;a { ) and the detector response distribution g 0 (v’;«i ) in equation (21 ) to indicate 
that this spectral index is inherited through the integral. 
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6. IDEAL DETECTOR 


The concept of a zero-resolution or ideal detector is very useful because it sets an upper bound 
on the expected performance of any real detector. Furthermore, it allows quantifing the magnitude of 
the uncertainty in the estimate of the spectral parameter. This uncertainty is measured in terms of the 
standard deviation of the estimator and attributable to event statistics (statistical fluctuation of inci- 
dent GCR proton energies) relative to the uncertainty in measuring the spectral parameter estimate 
induced by the detector’s nonzero energy resolution. 

Thus, for an ideal detector, p = 0 so that the standard deviation <r n£ = 0 for all GCR event 
energies E. Hence, the detector response to a GCR of energy E is given by Y = a + bE so that the 
incident energies may be directly obtained as E = (Y — a)/b so that the estimation procedures devel- 
oped in sections 4 and 5 apply. 

6.1 Method of Moments for a “Real” Detector 

The conditional expected value theorem is utilized that says the expected value of the condi- 
tional expected value is the unconditional expected value, 5 or in the notation of mathematical expec- 
tation applied to the detector response Y , 


/jy = <Y > = «Y \ E » . 

Thus, for a detector with constant resolution p, the following is obtained: 


( 22 ) 


jjy =(a + />p £ ) 


l+pp(p) j 
_ 1 / 


2 dx 


(23) 


where n Y is the mean of the detector response distribution and p £ is the mean of the simple power law 
distribution. The term involving the integral can be thought of as a correction term to the mean for the 
truncation and can be ignored whenever p < 0.30; i.e., 30-percent resolution or better. Using the 
method of moments, fd Y is estimated with the sample average Y . When combined with equation (5) 
with m = 1 for p £ yields equation (24) that can then be solved in terms of a, by numerical methods: 
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(24) 
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a | -1 




a\ 


- £^ _a| 


1 + 


pmp) J 


«i -2 


* v j 

e * or 


a \ _ a \ 


y[2n 


~/P 


For example, when the resolution is a constant 40 percent (p = 0.40), the point estimate of the 
spectral parameter a,, based on the 5,000 missions, is 2.801 using equation (24) and 2.79 using the 
same equation but with the correction term set to zero in the denominator. This results in a bias of 
= 0.01 that can be removed by including this correction term. This effect is much more pronounced 
when p = 0.6 and results in a bias of 0. 1 in the point estimate of a i , so the correction term is critical. 

When the detector response distribution is symmetric and truncation is negligible so that p Y 
= (a + bp E ), then a, can always be estimated using the mean of the detector responses Y to estimate 
p Y in equation (5) with m = 1 . This implies that knowledge of the variance of the detector distribution, 
and hence the resolution, is not required in order to estimate a t , provided it is known the resolution is 
<30 percent so the effect of truncation can be ignored. 

This is a useful result because if the uncertainty regarding the true resolution is non negligible, 
then the method of moments might be a good way to proceed with the estimation of or,; e.g., it is 
known the detector’s energy resolution is <30 percent but nothing more. However, as already noted, 
the method of moments does not provide the minimum variance estimator that the ML method does, 
which requires a complete specification of the detector parameters a, b, c, and d of this assumed 
Gaussian response function. Furthermore, the energy resolution of most real detectors is worse than 
30 percent. 

This estimator, based on the method of moments, is a function of the random variable Y and 
has its own associated pdf. Since Y has mean p Y and variance a Y , it is known by the central limit 
theorem the distribution of Y follows a normal distribution with mean p Y and variance rr Y /N. Thus, 
the variance of the detector response Y is (T^ = (t 2 ) - pp where 


( a ~ + 2 abp E + b 2 o E + b^ p E ^rj(p) 


7 0+p.r) 2 

\, V2? 


e 



(25) 


For example, when = 2.8, £, = 20 TeV, £ 2 = 5,500 TeV, and p = 0.40, Y is normally distributed 
with mean 131.58 GeV and standard deviation (213.69 GeV)/N^ 2 . The probability distribution of a j 
can be constructed and its mean and standard deviation obtained by solving the probability equation 
in equation (26) using the methods discussed with equation (11): 
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-J 


V2jc 


c 2 r/.v. (26) 


If the truncation effect is negligible in equation (25), then the following succinct formula for 
the variance of the detector response as a function of detector parameters a, b , and p and the mean p E 
and variance cr| of the power law distribution is obtained: 


0 



b 2 <7£ + p 2 (a + bp E y + b 2 o\ 


(27) 


In terms of the standard deviation of the detector response (7 V , the approximation in equation 
(27) is seen to be quite good. When p - 0.4, this formula yields o x = 213.37 GeV as compared to the 
exact value of 213.69 GeV obtained from equation (26) using the integral correction terms. When 
p = 0.6, this approximation yields G x = 237.31 GeV as compared to the actual value of 239.78 GeV. 
Thus, ignoring the truncation is not too serious when estimating the standard deviation but can be 
devastating for p > 0.4 when estimating the mean /Jy <ind hence oc^ when using the method of 
moments. Much insight into the estimation of the spectral parameter a, can be gleaned from equation 
(27) as it shows the relationship between the variance Gy of the detector response distribution, the 
variance g 2 e of the GCR proton energy spectrum, and the detector response function parameters a, b, 
and p. 


The influence of the variance and other higher moments of the simple power law energy spec- 
trum is visualized in figure 3, which shows the mean detector response (mean energy deposit) per 
mission for 30 simulated missions in comparison with the mean incident proton energy for 30 mis- 
sions. Corresponding standard deviations per mission are plotted in figure 4. Note that the detector 
response mean and standard deviation per mission tends to track the mean and standard deviation of 
the incident energies for the 30 missions, illustrating the strong influence of the GCR energy mission- 
to-mission fluctuations on the detector response variation, even in the presence of the “smearing” 
induced by this detector having 40-percent energy resolution. In section 6.2, the component of varia- 
tion due to the GCR event statistics will be the dominating the component of the total variation in the 
standard deviation of the estimator of the spectral index a,. 
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Figure 3. Comparing the mean incident energy with the mean 
of the detector responses for each of 30 missions. 
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Figure 4. Comparing the standard deviation of the GCR incident 

energies with the standard deviation of the detector responses 
for each of 30 missions. 


6.2 Maximum Likelihood for a “Real” Detector 

As in section 6. 1 , the method of ML seeks a ML which maximizes the log-L function so that log 
L(a ML )> log L(a,) for all a,, where the likelihood function for the detector response in the presence 
of the simple power law energy spectrum of N incident GCR protons over the energy range [E^E^ is 


logL(a, ) = 


N 

X lo s 


7=1 


N 

{$()(}' j' a \ ))J = X lo 8 

7=1 


E 2 

J }’{y j \E)(l) 0 (E:a l )clE 


(28) 


Because of the complexity of the integral and the desired capability to easily change the func- 
tional form of the detector response function g in equation (28), a numerical minimization algorithm 
called the Nelder-Mead downhill simplex method that does not require gradient information (deriva- 
tives) for obtaining a ML was chosen. 6 Since this is a minimization algorithm, the objective function is 
defined as 


N 

0(a,) = -logL(a 1 ) = -^log 

7=1 


E 2 

j g(yj I E)(f) 0 {E\a x )dE , 
E \ 


(29) 


so that minimizing 0(a ( ) maximizes log L(«j) as desired, where the integral is numerically evalu- 
ated. The following two termination criteria are used to halt the search procedure for the ML estimate 
at the (m + 1 )th iteration: 


(i) , ?!+ ] Kfj 

(ii) \0(a l m+l )-0(a Um )\<e 2 - (30) 

The search procedure continues until the termination criteria are met, which in words is as 
follows: (i) The movement in successive step sizes of oc ( is <£| and (ii) the objective function is 
changing by an amount <£ 2 . Typical values used for these two stopping tolerances are on the order of 
10 -5 and seem reasonable in light of the magnitude of the parameter being estimated (=2.8) and the 
value of the objective function in the vicinity of the found ML solution 0(ct ^ L ), of the order of 
magnitude 1 0 5 when E | is taken to be anywhere between 1 0 and 30 TeV, so the number of terms in the 
sum is between 182,000 to 26,000, respectively. Furthermore, changing £, and/or e 2 in either direc- 
tion by an order of magnitude provided no noticeable change in results. 
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Figure 5 shows the ML estimates of a, for a zero-percent resolution detector obtained from 
equation (14) in comparison with the ML estimates obtained from a 40-percent resolution detector 
and applying the downhill simplex algorithm to equation (29) for 30 missions. This very close com- 
parison suggests that the GCR event statistics are the dominating component of uncertainty in the 
estimation of the spectral parameter 



Mission 

Figure 5. Maximum likelihood estimates for zero- and 40-percent 
resolution detector for 30 missions. 


18 


7. NON-GAUSSIAN DETECTOR RESPONSE FUNCTION 


The results so far have assumed a Gaussian detector response function. While reference 4 
suggests that a Gaussian function is reasonable, there is concern that perhaps the response function is 
skewed slightly to the right and that this “tail” will contribute to greater difficulties in estimating the 
spectral parameter. The gamma response function, capable of describing a wide variety of shapes with 
right-hand skewness (outer curve from the right in fig. 6), was introduced to address this concern. Its 
parameters were set to provide a constant 40-percent energy resolution over the study range 20-5,500 
TeV from which events having incident energy described by a simple power law with a x = 2.8 were 
simulated. The number of these events was determined using the baseline detector collecting power. 



Energy Deposit (GeV) 

Figure 6. Gaussian and gamma detector response functions 
to 40 TeV proton (40-percent resolution). 


Based on 10,000 simulated missions, the mean ML estimate of a, was 2.800 and standard 
deviation 0.0093 as compared to 2.800 and 0.0092 for the mean and standard deviation, respectively, 
corresponding to a Gaussian response function. Thus, the conclusion that this skewness and other 
non-Gaussian assumptions are insignificant when using the method of ML. However, as will be shown, 
the real key is not so much a matter of what the response function is as knowing what it is; i.e., an 
accurately understanding of the response function over the incident energy range must be made. It 
should also be noted that while the gamma response function as modeled has a constant energy reso- 
lution of 40 percent, the Gaussian used in this comparison had a 39-percent resolution as previously 
discussed. 
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8. ENERGY-DEPENDENT RESOLUTION STUDY 


The situation in which the detector response function is assumed to be Gaussian, but its energy 
resolution varies with incident GCR event energy, is of particular interest to designers of cosmic-ray 
detectors. In the studies presented so far in this TP, the detector response function is assumed to be 
Gaussian with a linear mean response (energy deposit) of the form ( a + bE) and with constant detector 
energy resolution p so that the parameter c in the Gaussian response function is defined as a(E) 
=p{a + bE). Two cases of interest are as follows: ( I ) Energy resolution is “getting better,” from 
40-percent resolution at E i =20 TeV to 30 percent at E 2 = 5,500 TeV and (2) “getting worse,” from 
30-percent resolution at E i =20 TeV to 40 percent at E 2 = 5,500 TeV. These two cases are modeled by 
assuming that o{E) is a linear function of incident GCR energy of the form (c + dE) and then the 
coefficients c and d are determined by matching the conditions for each of the two cases. Doing so 
yields the energy-dependent resolution curves depicted in figure 7. 

Table 1 shows the results based on 100 simulated missions using the same incident GCR ener- 
gies for both cases. The mean estimates are seen as essentially unbiased with standard deviations 
having expected comparisons; e.g., standard deviations slightly larger for the “getting worse” case. 
The constant 30- and 40-percent cases are included for comparison. 


Nonconstant Resolution Curves 
tor Energies Between 20-5,500 TeV 



Figure 7. Energy-dependent resolution curves. 
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Table 1 . Nonconstant energy resolution study. 


Mean and Standard Deviation of the Maximum Likelihood Estimate ol ai 
Based on 5,000 Missions 

Spectral 

Parameter 

Resolution 

Constant 30% 

Constant 40% 

Nonconstant 
(Getting Better) 

Nonconstant 
(Getting Worse) 

ai 

Mean Std. Dev. 
2.8 0.0088 

Mean Std. Dev. 
2.8 0.0092 

Mean Std. Dev. 
2.8 0.0090 

Mean Std. Dev. 
2.8 0.0091 




9. IMPLICATIONS OF DETECTOR RESPONSE MODEL UNCERTAINTIES 


Maximum likelihood estimation of the spectral parameter a 1 using equation (29) requires the 
complete specificity of all detector response model parameters. The reality of knowing these param- 
eters with little or no surrounding uncertainty depends largely on designers being able to calibrate the 
detector at different incident energies at a particle accelerator facility. However, because space-based 
detectors will be exposed to GCR events having energy much greater than those energies available at 
accelerator facilities, it becomes essential to gain an understanding of the detector's response func- 
tion using Monte Carlo simulations of the detector’s response (energy deposit) to those energies that 
cannot be attained at accelerator facilities. These simulations, coupled with a favorable comparison 
between simulation results and accelerator results at energies available in a test facility, will provide 
a better understanding of the detector response function. 

Next, by way of example, an investigation is done of the impact on estimating a | when certain 
detector response function parameters are incorrectly known, and as shall be seen, this state of igno- 
rance will manifest itself as a bias in the mean or point estimate of a f . This state of ignorance is 
modeled by simulating detector responses according to one set of detector response function param- 
eters and then using a different set of parameters in the detector response function g in equation (29) 
during the ML estimation procedure. 

Since detector resolution is an important design parameter, consider first the case where the 
detector has a constant energy resolution. In an assumed state of misunderstanding, a different resolu- 
tion value was used in equation (29). For example, suppose the real detector resolution is a constant 
35 percent, but in the simplex search, the resolution parameter p is set to different constant values in 
equation (29), corresponding to resolutions ranging from 30 to 40 percent. This situation is modeled 
by simulating the detector responses K f - as 

Yj = {a + bEj)(\ + 0.35 Z f ) (31) 

according to equation (17) and then set p to the different values in equation (29) in the ML procedure. 
The results for 500 simulated missions for each assumed resolution 30 through 40 percent while the 
“real” resolution is a constant 35 percent are presented in table 2. Note that the mean estimates exhibit 
a bias as a result of using incorrect values of p in equation (29). Note when p = 0.35 in equation (29) 
and hence matches the “real” resolution as used in equation (31) to simulate the detector responses, 
the means of the ML estimates match the assumed spectral parameter value of a ] = 2.8 used in the 
simulation. 
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Table 2. Biased estimate of a, when detector energy resolution is incorrectly 
known and assumed to be a constant 35 percent. 


Real Energy 
Resolution 

30% 31% 32% 33% 34% 35% 36% 37% 38% 39% 40% 

Mean Estimate 

Of Oli 

2.83 2.83 2.82 2.81 2.81 2.80 2.79 2.79 2.78 2.77 2.76 


Another interesting situation to consider is when the real detector resolution is energy depen- 
dent, but in a state of ignorance, a constant resolution value is used in equation (29) in the ML estima- 
tion of CC|. For example, if the real detector resolution is “getting better as in figure 7 but instead a 
constant p = 0.40 is used in equation (29), a mean estimate of 2.81 for a.\ is obtained, whereas it p = 
0.35 is used in equation (29), a mean estimate of 2.78 is obtained. One may initially be surprised by 
this comparison since 0.35 can roughly be considered the average resolution value for this energy- 
dependent case and that using a constant p = 0.35 in equation (29) should provide a closer result (less 
bias) than when p is set to 0.40. However, remember that the p = 0.40 assumption is closest to “real- 
ity” in the region where events are most numerous because of the steepness of the power law and the 
true energy-dependent case does start at 40-percent resolution. Also interesting is that if p= 0.38 in 
equation (29), the mean estimate 2.8 is obtained, even though the “real resolution is the energy- 
dependent “getting better” case— a result discovered by trial and error! 

A much larger bias in the estimate of Cfj is observed when the case where the real detector 
response function is the gamma function with a constant energy resolution ot 40 percent, as depicted 
in figure 6, but instead, in a state of misunderstanding, a Gaussian response function in the ML proce- 
dure is assumed and the Gaussian response function defined in equation (18) is substituted into equa- 
tion (29), resulting in a mean estimate of 2.52 for a,. 
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10. APPLICATION TO REAL COSMIC-RAY DATA 


In these simulations, the GCR events are simulated from an energy range £ j to E->, where 
typically E 2 = 5,500 TeV for this generic-sized detector and £, is a value between 5 and 25 TeV. The 
choice of£ 2 is based on the collecting power of the detector and is chosen such that there will be only 
a negligible number of events above £ 2 . The selection of £ ( is largely dictated by the practical number 
of events that can be handled in the simulation for several thousand missions. 

Next, for each of these simulated GCR events, a detector response is simulated according to 
the assumed detector response function and then the full set of simulated responses are used to esti- 
mate the spectral parameters. However, because no energies below E l are simulated, frequency histo- 
grams of the simulated detector responses do not match the front-end portion of a real cosmic-ray 
energy spectrum. This difference or mismatch is an artifact of not having generated events from 
below £, that would have had the effect of filling in this front-end portion of the histogram and for 
consequently resembling a real cosmic-ray energy spectrum. 

This difference is not critical when making relative comparisons of the effects of design 
parameters when the detector response function parameters used to generated the simulated responses 
match those detector response function parameters used in equation (29) in the simplex search for the 
ML estimate of a i.e., a perfect understanding is had of the detector response function. However, 
when the impacts of response function uncertainties are investigated, it is more important that the 
simulation techniques produce results that are closer to a real cosmic-ray energy response spectrum. 

To make the histogram of the simulated detector responses look like a real cosmic-ray energy 
spectrum, a cut v ( . in the simulated detector responses is introduced and all energy deposits smaller 
than y c are dropped. In the simulation, the choice of y c dictates admissible values of E i because E > 
must be chosen so that only a negligible number of events having incident energy <£, deposit ener- 
gies >y t ., which of course depends on the detector’s energy resolution. For example, if y = 80 GeV 
corresponding to the mean energy deposit of a 25-TeV GCR event and a Gaussian response function 
having 40-percent energy resolution and mean response ( a + bE ) is considered, as used for the baseline 
detector and defined in equation (18), then £] can be any value <10 TeV, since only a negligible 
number of events from below' 10 TeV will deposit more than 60 GeV. Selecting £j = 10 TeV provides 
~ 1 82,000 GCR events and setting y ( . = 80 GeV and dropping all simulated detector responses smaller 
than y c leaves =44,500 events on average above y c , which produces a simulated response spectrum 
that does indeed resemble a real response spectrum. Table 3 shows the average number of events with 
energy deposits >y { . for selected values of £, from 2 to 20 TeV. The practical benefit of table 3 is that 
it provides the optimal value to set £ j in the simulation for a given cut y so as to minimize the number 
of events to initially simulate and achieve the same result. For example, if the cut is set at 60 GeV. 
then £| =4 TeV would be set, since there is no additional benefit by taking £| < 4 TeV. 
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Table 3. Number of simulated energy deposits above cut y ( . for various E , 

for Gaussian detector response function with 40-percent resolution. 



Cut y c 
(GeV) 

54 

59 

67 

75 

80 

93 

106 

«W>£i) 
3,293,737 
1,587,539 
945,876 
632,988 
455,899 | 

345,433 
271,630 
219,737 
181,777 
153,119 
130,921 
113,354 
99,198 
87,613 
78,004 
69,939 
63,101 
57,249 
52,200 

*1 (TeV) 
2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

127,707 

127,605 

126,796 

124,971 

121,685 

116,887 

110,969 

104,270 

97,235 

90,162 

83,311 

76,840 

70,819 

65,292 

60,242 

55,655 

51,465 

47,714 

44,320 

97,784 

97,828 

97,648 

97,196 

96,151 

94,211 

91,438 

87,799 

83,592 

78,980 

74,209 

69,427 

64,801 

60,377 

56,217 

52,371 

48,748 

45,464 

42,451 

69,979 

70,013 

69,936 

69,901 

69,731 

69,358 

68,644 

67,457 

65,793 

63,680 

61,201 

58,464 

55,625 

52,713 

49,817 

47,008 

44,292 

41,707 

39,295 

52,896 

52,906 

52,846 

52,841 

52,820 

52,747 

52,621 

52,313 

51,766 

50,982 

49,866 

48,500 

46,941 

45,202 

43,357 

41,454 

39,562 

37,678 

35,835 

44,873 

44,878 

44,837 

44,812 

44,808 

44,779 

44,769 

44,646 

44,412 

44,073 

43,463 

42,658 

41,688 

40,518 

39,237 

37,823 

36,396 

34,923 

33,436 

31,487 
31,487 
31,468 
31,451 
31 ,447 
31,426 
31,461 
31,440 
31,416 
31,404 
31,313 
31,146 
30,907 
30,576 
30,138 
29,583 
28,967 
28,276 
27,510 

23,440 

23,465 

23,434 

23,424 

23,414 

23,407 

23,427 

23,407 

23,389 

23,427 

23.418 

23,381 

23,352 

23,286 

23,165 

23,003 

22,785 

22,528 

22,206 


A very important benefit is realized by introducing the cut y c ; that is, the lower limit of integra- 
tion in equation (29) can be set to any value E ^ < E j with the implication that this ML procedure can 
be made independent of the range of integration, as long as E L is chosen wisely using table 3 (or its 
equivalent for detectors having a different energy resolution), so that the ML estimation procedure 
herein developed can now be applied to real cosmic-ray detector response data. It should be noted that 
cuts on the high end are not necessary to consider since any value E ^ > E 2 is suitable because the 
number of events above £-> is negligible and so </>q is approximately zero when E > £3 the integral 
definition of the detector response pdf. However, setting E H unnecessarily high would result in many 
unnecessary calculations in the numerical integration of equation (29). 


Introduction of a cut y\. requires a modification to the objective function in equation (29) to 
handle the conditional detector response distribution because of the constraint y > _y r Thus, the objec- 
tive function using cut y c becomes 


0(a,)=- 



(32) 
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where 


\^" g(y, I E:p)(f> 0 (E\a { )dE 

^(>-1 v- >.v ( . ;«,) = — — . >'/ > >V • (33) 

i - J £o(v; «i )dy 

o 

From a simulation point of view, E j = 2 TeV is about the lowest value that was used because of 
the vast number of generated events and the requirement to handle thousands of simulated missions 
necessary for meaningful inferences. Consequently, cuts much less than 60 GeV were not used in 
these studies. However, cuts in real cosmic-ray data can be taken to be much lower since the real 
spectrum is automatically Filled in from events having incident energies much less than 2 TeV and the 
lower limit of integration can then be taken to be any small value. 
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11. CONFIDENCE INTERVAL FOR a x 


The ML estimate a ML obtained from equation (29) was shown to (1) be unbiased, (2) attain the 
Cramer-Rao minimum variance bound, and (3) be normally distributed when the number ot events N is 
= 1,200 or larger. Thus, the variance of a ML is Var(a ML ) = -1 l[d 2 LLIda\] evaluated at a ML for any given 
“mission” and where LL denotes the log-likelihood function given in equation (28) which can readily be 
obtained by numerical methods. A confidence interval can then be constructed for the unknown spectral 
parameter a as illustrated for each of 20 missions at 95-percent confidence in figure 8. 

Extending the number of simulated missions to 2,000 provides a frequency histogram of a ML as 
shown in figure 9 and of the estimated standard deviation of O^l as shown * n figure 10. 


95-Percent Confidence Interval for a 

— ♦— Lower 95-Percent Confidence Limit — ■— a ML — ±— Upper 95-Percent Confidence Limit 



Figure 8. Maximum likelihood estimate with 95-percent confidence interval 
for 20 “missions.” 
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Frequency 




Figure 10. Frequency histogram of the estimated standard 
deviation of a ML (2,000 missions). 



12. TESTING FOR SLOPE DIFFERENCES OF TWO COSMIC-RAY ELEMENTAL SPECIES 


Given two cosmic-ray elemental species A and B with slopes oc and (3, an important hypothesis to 

test is 


Hq. a- /3 = 0 (same “slopes”) 

versus 

//,: a-13* 0 . 

To test this hypothesis, ML estimates a ML and (3^ for each slope parameter using equation (29) and their 
respective variances Varfc:^) an d ^ ar fAvlL^ numer i ca ^y obtained as discussed in section 1 1, Because 
a ML and are both normally distributed for large N, the test statistic T = (a ML - is used to test H 0 . 

Note the following: 

1 . T has mean fi T = (a- /?) because a ML ar| d are unbiased. 

2. The standard deviation of T is (J T = + Var (/3 M , ) because a ML and ^ are independent. 

3. T is normally distributed because a ML and /3 ML follow the normal distribution. 

4. T is normal(0, Oj) when Hq is true. 

A two-sided test rejects Hq with confidence C whenever ITI > w h ere -c/2 * s cr ^' ca ^ va l ue 

of the normal distribution for confidence level C for a two-sided test. The width of the confidence interval 

for (a - P) is then 2 

Furthermore, if the number of events is doubled, then let D denote the test statistic similarly defined 
as T. Then D is normally distributed with mean zero and standard deviation <J D under H 0 also; note that C D 
= o T /Sqrt(2) since the variance of each ML estimate is halved when the number of events doubles as 
discussed in section 13. Thus, the width of the corresponding confidence interval is reduced by a factor of 
1/Sqrt(2) when the number of events doubles. 
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13. SUMMARY, REMARKS, AND CONCLUSIONS 


Two methods for estimating the single spectral index a t of a simple power law have been 
investigated. The first method, called the method of moments, was found to be very useful in studying 
the general nature of the statistical estimation problem as well as yielding an analytical solution that 
could be compared with Monte Carlo simulation results. Furthermore, when the detector resolution is 
better than 30 percent so that the truncation of the detector response function is negligible, the method 
of moments provides an estimator of a j without requiring specific knowledge of the detector resolu- 
tion p but only that it is better than 30 percent. This does not imply p is insignificant when it is 
<30 percent, but only that the correction terms in equations (24)-(26) can be ignored. Thus, explicit 
knowledge of the value of p is not needed to estimate oq. In fact, the standard deviation of the estima- 
tor increases as p increases, as one would expect, and results from the fact that whatever phappens to 
be, its impact is communicated to the estimate of through the variance of the detector mean 
response Y — a function of p as indicated in equations (26)— (28). Another interesting result is that 
when the resolution is <30 percent, the explicit functional form of the detector model need not be 
known, but only that it is symmetric. Unfortunately, most detector response functions are worse than 
30-percent resolution and may be asymmetric as well. 

The method of ML estimation clearly stands out as the method of choice for estimating a j in 
terms of minimum variance and consistency (asymptotically unbiased) as well as asymptotic normal- 
ity that allows for probabilistic statements, such as confidence intervals for the unknown spectral 
parameter. These results as a function of detector resolution are shown in figure 8. 

Maximum Likelihood and Method of Moments 
Estimator of Spectral Parameter a, 

(Simple Power Law) Versus Detector Resolution 


2. Standard Deviation, Method of Moments (Simulation) 
— B— 4. Standard Deviation, Maximum Likelihood (Simulation) 



Figure 1 1 . Comparison between method of moments 
and ML as a function of detector resolution. 
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When compared to the standard deviation of the method of moments estimator, the ratio varies 
from 1 .47 for the zero-percent resolution detector to 1 .33 for the 50-percent resolution detector, which 
is roughly equivalent to losing half of the detector’s collecting power by choosing the interior method 
of moments estimation technique. 

Also shown is that the standard deviation of the estimate for both estimation procedures is 
inversely proportional to the square root of the sample size, so that halving the collecting power 
increases the standard deviation by a factor of 4l . This holds true for the standard deviation of the 
ML estimate as long as it attains the Cramer-Rao lower bound, which it does when the number o 

GCR events exceeds 2,000. 

Another important result is the relationship between the collecting power and the energy reso- 
lution of the detector. A measure of the detector’s ability to estimate the spectral parameter a, is its 
standard deviation, as seen in figures 8 and 9. The dominant component of the standard deviation ot 
~ is attributable directly to the large fluctuations in GCR incident energies, being driven by the 
large variance and other higher moments of the simple power law distribution. This large component 
can only be reduced by increasing the number of events N that is controlled by the collecting power ot 
the detector. A comparison of the standard deviation of a ML for the generic detector discussed in this 
TP and when its collecting power is halved is shown in figure 9. Table 4 provides the numerical 
results used to construct many of the figures in this section. 


Collecting Power Comparison: Halt as Many Events 


Maximum Likelihood (/V=52.000) 
Maximum Likelihood (/V=26,000) 



Figure 12. Comparing the effect of collecting power on the standard 
deviation of the ML estimate of the spectral index a, . 
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Table 4. Numerical values used to construct figures 1 1 and 12 GCR. 


£|=20 TeV, E 2 =5,500 TeV, 

Detector Resolution 





<*1=2.8, ^average=52,000 Events. 

5,000 Mission Averages tor Simulation Results. 

0% 

20% 

40% 

50% 

1 . Method of moments (theory) 

0.0115 

0.0116 

0.0128 

0.0136 

2. Method of moments (simulation) 

0.0114 

0.0117 

0.0125 

0.0133 

3. Maximum likelihood (Cramer-Rao lower bound) 

0.00786 

Analytica 

solution not 

available 

4. Maximum likelihood (simulation) 

0.0078 

0.0083 

0.0092 

0.0100 

5. Mean detector response (GeV) (theory) 

130.66 

130.66 

131.58 

138.85 

6. Mean detector response (GeV) (simulation) 

130.66 

130.64 

130.64 

138.81 

7. Standard deviation (theory) 

192.07 

197.61 

213.69 

239.77 

8, Standard deviation (simulation) 

191.47 

196.86 

213.33 

238.82 

9. Coefficient of variation Vy (detector, %) 

147 

151 

162 

173 

Ei =20 TeV, E 2 =5,500 TeV, o 1 =2.8, « awaBe =26,000 
Events. 5,000 Mission Averages for Simulation Results. 





10. Maximum likelihood 

0.0110 

0.0118 

0.0132 

0.0144 

1 1 . Ratio of line 4 to line 1 0. compare to V ~ 

1.41 

1.42 

1.43 

1.44 
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